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1 Introduction 

On the nanoscale, materials around us have a surprisingly simple structure: The standard model 
of solid state physics and chemistry only knows of two types of particles, namely the nuclei 
making up the periodic table and the electrons. Only one kind of interaction between them needs 
to be considered, namely the electrostatic interaction. Even magnetic forces are important only 
in rare occasions. All other fundamental particles and interactions are irrelevant for chemistry. 
The behavior of these particles can be described by the Schrodinger equation (or better the 
relativistic Dirac equation), which is easily written down. However, the attempt to solve this 
equation for any system of interest fails miserably due to what Walter Kohn termed the expo- 
nential wall [[II. 

To obtain an impression of the powers of the exponential wall, imagine the wave function of a N2 
molecule, having two nuclei and fourteen electrons. For particles, the Schrodinger equation 
is a partial differential equation in 3A^ dimensions. Let us express the wave function on a grid 
with about 100 points along each spatial direction and let us consider two spin states for each 
electron. Such a wave function is represented by 2^'*100^*^^ ^ 10^°° complex numbers. A data 
server for this amount of data, made of current terabyte hard disks, would occupy a volume with 
a diameter of 10^° light years! 

Treating the nuclei as classical particles turned out to be a good approximation, but the quantum 
nature of the electrons cannot be ignored. A great simplification is to describe electrons as non- 
interacting quasi particles. Instead of one wave function in 3A^ dimensions, one only needs to 
describe N wave functions in three dimensions each, a dramatic simplification from 10^°° to 
10^ numbers. 

While the independent-particle model is very intuitive, and while it forms the basis of most 
text books on solid-state physics, materials physics, and chemistry, the Coulomb interaction 
between electrons is clearly not negligible. 

Here, density-functional theory [[2l[3][ comes to our rescue: it provides a rigorous mapping from 
interacting electrons onto a system of non-interacting electrons. Unfortunately, the exact map- 
ping is utterly complicated and this is where all the complexity goes. Luckily, there are simple 
approximations that are both, intuitive and surprisingly accurate. Furthermore, with the help of 
clever algorithms, density-functional calculations can be performed on current computers for 
large systems with several hundred atoms in a unit cell or a molecule. The microscopic insight 
gained from density functional calculations is a major source of progress in solid state physics, 
chemistry, material science, and biology. 

In the first part of this article, I will try to familiarize the novice reader with the basics of density- 
functional theory, provide some guidance into common approximations and give an idea of the 
type of problems that can be studied with density functional theory. 

Beyond this article, I recommend the insightful review articles on density functional theory 
by Jones and Gunnarsson [[3, Baerends [[3, von Barth [[6|[, Perdew [[7][, Yang [[H, and their 
collaborators. 

Solving the one-particle Schrodinger equation, which results from density-functional theory. 
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for real materials is a considerable challenge. Several avenues have been developed to their 
solution. This is the field of electronic structure methods, which will be discussed in the second 
part of this article. This part is taken from earlier versions by Clemens Forst, Johannes Kastner 
and myself [HI [la. 

2 Basics of density-functional theory 

The dynamics of the electron wave function is governed by the Schrodinger equation ihdt \^) = 
H]"^) with the N-particle Hamiltonian H. 



With rrie we denote the electron mass, with cq the vacuum permittivity, e is the elementary 
charge and h is the Planck quantum divided by 27r. The Coulomb potentials of the nuclei have 
been combined into an external potential Vextir). 

All N-electron wave functions ^{xi, . . . ,xn) obey the Pauli principle, that is they change their 
sign, when two of its particle coordinates are exchanged. 

We use a notation that combines the position vector f E M'^ of an electron with its discrete 
spin coordinate a E {t, i} into a single vector x := (r, o"). Similarly, we use the notation of 
a four-dimensional integral J d^x := J2a I ^^"^ ^^^^ ^P^'^ indices and the integral 

over the position. With the generalized symbol 5{x — x') := 5ay5{r — r') we denote the product 
of Kronecker delta of the spin coordinates and Dirac's delta function for the positions. While, 
at first sight, it seems awkward to combine continuous and discrete numbers, this notation is 
less error prone than the notation that treats the spin coordinates as indices, where they can be 
confused with quantum numbers. During the first reading, the novice will ignore the complexity 
of the spin coordinates, treating x like a coordinate. During careful study, he will nevertheless 
have the complete and concise expressions. 

One-particle reduced density matrix and two-particle density 

In order to obtain the ground state energy E = {^IH]^) we need to perform 2^ integrations in 
3A^ dimensions each, i.e. 



However, only two different types of integrals occur in the expression for the energy, so that 
most of these integrations can be performed beforehand leading to two quantities of physical 
significance. 

• One of these quantities is the one-particle reduced density matrix p^^^ (x, x') , which allows 
one to evaluate all expectation values of one-particle operators such as the kinetic energy 




(1) 



E 
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and the external potential, 



The other one is the two-particle density n^'^\r, r'), which allows to determine the inter- 
action between the electrons, 

n^^\r,r') ■= N{N - 1) } ^ I d^X3 . . . I d^XN \'^{x,x',X3, . . . ,xn)\'' ■ (4) 



N{N - 1) / ^^^3 ■ ■ ■ j d^XN 1^ 



If it is confusing that there are two different quantities depending on two particle coordinates, 
note that the one-particle reduced density matrix p^^^ depends on two x-arguments of the same 
particle, while the two-particle density n^'^^ depends on the positions of two different particles. 
With these quantities the total energy is 

E = J d^x' j d^x 5{x' - x) {^^^ + Ve^t{r)^ P^^^ (x, x') 

+ - / d'r / d\' , (5) 

^ J J 47reo|r — r'\ 

where the gradient of the kinetic energy operates on the first argument r of the density matrix. 

One-particle reduced density matrix and natural orbitals 

In order to make oneself familiar with the one-particle reduced density matrix, it is convenient 
to diagonalize it. The eigenstates ipn{r) are called natural orbitals [fTT]| and the eigenvalues 
are their occupations. The index n labels the natural orbitals may stand for a set of quantum 
numbers. 

The density matrix can be written in the form 

= Y,fn^n{S)^n{^') . (6) 

n 

The natural orbitals are orthonormal one-particle orbitals, i.e. 

d'^X Lp*^{x)(fn{x) = 5m,n ■ (V) 

Due to the Pauli principle, occupations are non-negative and never larger than one lfT2l . The 
natural orbitals already point the way to the world of effectively non-interacting electrons. 
The one-particle density matrix provides us with the electron density 

n«(f) = 5^p(in^,^) = EE/"^n(^)¥'n(^) . (8) 

a an 

With the natural orbitals, the total energy Eq.[5]obtains the form 

E = Y,fnj d^x )^ VVn(x) + j dh i;ext(r)n(^)(f) 

+ - dh d^r' . (9) 

2 J J 47reo|r — r'l 
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Two-particle density and exchange-correlation hole 

The physical meaning of the two-particle density n^'^\r,r') is the following: For particles 
that are completely uncorrected, meaning that they do not even experience the Pauli princi- 
ple, the two particle density would be0 the product of one-particle densities, i.e. n^^\r,r') = 
n^^\r)n^^\r') . If one particle is at position tq, the density of the remaining — 1 particles is 
the conditional density 

n«(fo) • 

The conditional density is the electron density seen by one of the electrons at tq. This observer 

electron obviously only sees the remaining — 1 electrons. 

It is convenient to express the two-particle density by the hole function h{f, r'), i.e. 



n^^\r') + h{r,r') 



(10) 



One electron at position r does not "see" the total electron density n^^^ with N electrons, but 
only the density of the — 1 other electrons, because it does not see itself. The hole function 
h{r\j, r) is simply the difference of the total electron density and the electron density seen by 
the observer electron at tq. 

The division of the two-particle density in Eg. [TOl suggests to split the electron-electron interac- 
tion into the so-called Hartree energy 

Eh = - I <rr I d'\ — ^-^ (11) 



2 J J 47reo|r 
and the potential energy of exchange and correlation 

J J 47reo|r — r'\ 

Keep in mind that Uxc is not the exchange correlation energy. The difference is a kinetic energy 
correction that will be discussed later in Eq.[T9l 

The hole function has a physical meaning: An electron sees the total density minus the electrons 
accounted for by the hole. Thus each electron does not only experience the electrostatic poten- 
tial of the total electron density n^^\r), but also the attractive potential of its own exchange 
correlation hole h{fo, r). 
A few facts for this hole density are apparent: 

1. Because each electron of a N-electron system sees N— 1 other electrons, the hole function 
integrates to exactly minus one electron 

£rh{fQ,r) = -l (13) 

irrespective of the position tq of the observing electron. 



' This is correct only up to a term that vanishes in the limit of infinite particle number. 
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Fig. 1: Exchange hole in silicon. The cross indicates the position of the observer electron. The 
black spheres and the lines indicate the atomic positions and bonds in the (110) plane. Reprinted 
figure with permission from Mark S. Hybertsen and Steven G. Louie, Physical Review B 34, 
5390 (1986). Copyright 1986 by the American Physical Society 



2. The density of the remaining — 1 electrons can not be larger than the total electron 
density. This implies 

h{fo,f) > -n''^\r} . (14) 

3. Due to the Pauli principle, no other electron with the same spin as the observer electron 
can be at the position tq. Thus the on-top hole h{fQ, tq) obeys the limits [[T3l 

- In^'Hro) > h{fo,fo) > -n«(fo) . (15) 

4. Assuming locality, the hole function vanishes at large distances from the observer electron 
at ro, i.e. 

h{fo,r) — for |r — fo\ ^ oo . (16) 

With locality I mean that the density does not depend on the position or the presence of 
an observer electron, if the latter is very far away. 



A selfmade functional 

It is fairly simple to make our own density functional^: For a given density, we choose a simple 
shape for the hole function, such as a spherical box. Then we scale the value and the radius such 
that the hole function integrates to —1, and that its value is opposite equal to the spin density at 
its center. The electrostatic potential of this hole density at its center is the exchange-correlation 
energy for the observer electron. Our model has an exchange correlation energ)§ of 

functional F[y\ maps a function y{x) to a number F. It is a generalization of the function F{y) of a vector 
y, where the vector index of y is turned into a continuous argument x. 

■^For this model we do not distinguish between the energy of exchange and correlation and its potential energy 
contribution 
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Fig. 2: Left: Scheme to demonstrate the construction of the exchange correlation energy from 
a simple model. Right: exchange correlation energy per electron e^c as function of electron 
density from our model, Hartree-Fock approximation and the exact result. The symbol "Na" 
indicates the density of Sodium. 

The derivation is an elementary exercise and is given in the appendix. The resulting energy 
per electron e^c is given on the right-hand side of Fig. [2] indicated as "model" and compared 
with the exact result indicated as "LSD" and the Hartree-Fock result indicated as "HF" for a 
homogeneous electron gas. 

The agreement with the correct result, which is surprisingly good for such a crude model, pro- 
vides an idea of how robust the density-functional theory is with respect to approximations. 
While this model has been stripped to the bones, it demonstrates the way physical insight enters 
the construction of density functionals. Modem density functionals are far more sophisticated 
and exploit much more information [[T4|. but the basic method of construction is similar. 



Kinetic energy 

While the expression for the kinetic energy in Eq. |9] seems familiar, there is a catch to it. In 
order to know the natural orbitals and the occupations we need access to the many-particle 
wave function or at least to its reduced density matrix. 

A good approximation for the kinetic energy of the interacting electrons is the kinetic energy 
functional Ts[n^^^ of the ground state of non-interacting electrons with the same density as the 
true system. It is defined by 



Ts[n 



(1)1 



mm 

{/ne[0,l],|»/'n>} 



+ / d'r Veffif^ 



2m 



XI 



(18) 



Note that /„ 7^ /„ and that the so-called Kohn-Sham orbitals ipni^) diffei0 from the natural 
''To be precise, Kohn-Sham orbitals are the natural orbitals for non-interacting electrons of a given density. 
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orbitals fn{x). Natural orbitals and Kohn-Sham wave functions are fairly similar, while the 
occupations /„ of Kohn-Sham orbitals differ considerably from those of the natural orbitals. 
The effective potential Veff{r) is the Lagrange multiplier for the density constraint. „ is the 
Lagrange multiplier for the orthonormality. Diagonalization of A yields a diagonal matrix with 
the one-particle energies on the diagonal. 

This kinetic energy T5[n'^^)] is a unique functional of the density, which is the first sign that we 
are approaching a density-functional theory. Also it is the introduction of this kinetic energy, 
where we made for the first time a reference to a ground state. Density functional theory as 
described here is inherently a ground-state theory. 

Why does the true kinetic energy of the interacting system differ from that of the non-interacting 
energy? Consider the hole function of a non-interacting electron gas. When inserted into Eq.[T2l 
for the potential energy of exchange and correlation, we obtain a contribution to the total 
energy that is called exchange energy. The interaction leads to a second energy contribution 
that is called correlation energy. Namely, when the interaction is switched on, the wave func- 
tion is deformed in such a way that the Coulomb repulsion between the electrons is reduced. 
This makes the hole function more compact. However, there is a price to pay when the wave 
functions adjust to reduce the Coulomb repulsion between the electrons, namely an increase of 
the kinetic energy: Pushing electrons away from the neighborhood of the reference electrons 
requires to perform work against the kinetic pressure of the electron gas, which raises the ki- 
netic energy. Thus, the system has to find a compromise between minimizing the electrostatic 
repulsion of the electrons and increasing its kinetic energy. As a result, the correlation energy 
has a potential-energy contribution and a kinetic-energy contribution. 

This tradeoff can be observed in Fig. [21 The correct exchange correlation energy is close to our 
model at low densities, while it becomes closer to the Hartree-Fock result at high densities. This 
is consistent with the fact that the electron gas can easily be deformed at low densities, while 
the deformation becomes increasingly costly at high densities due to the larger pressure of the 
electron gas. 

The difference between and the true kinetic energy is combined with the potential energy of 
exchange and correlation Uxc from Eq. [121 into the exchange correlation energy E^c, i-C- 



Note, that the (pni^) and the /„ are natural orbitals and occupations of the interacting electron 
gas, and that they differ from the Kohn-Sham orbitals ^n{x) and occupations /„. 



They are however different from the natural orbitals of interacting electrons at the same density. 




n 



(19) 
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Total energy 

The total energy obtains the form 



E = min [S" f d^x ^l{x)-^V^i)n{S) 



+ / d-^r Veff{r) 



x] 



n{f) ) + / di^r Vextif)n^^\'r) 



+ \ I d'r ! dh' ^'^^'HyO^^^HO ^ _ V- ^^^^ ('(^„|^„,) - ^^^^'j 1 . (20) 
2i J 47reo|r-r'| TT;^ V /J 



In order to evaluate the total energy with Eq. |20l we still have to start from the many-particle 
wave function |$). Only the many -particle wave function allows us to evaluate the one-particle 
density n^^\'r) and the exchange correlation energy E^c- Kohn-Sham orbitals and occupa- 
tions fn are obtained by an independent minimization for each density. 

If, however, we were able to express the exchange-correlation energy E^^. as a functional of the 
density alone, there would be no need for the many-particle wave function at all and the terrors 
of the exponential wall would be banned. We could minimize Eq.|20]with respect to the density, 
Kohn-Sham orbitals and their occupations. 

Let us, for the time being, simply assume that E^dn^^^] is a functional of the electron density 
and explore the consequences of this assumption. Later, I will show that this assumption is 
actually valid. 

The minimization in Eq. [20] with respect to the one-particle wave functions yields the Kohn- 
Sham equations 



+ Veffif) - e. 



ll)n{x) = with / d'^X 1prn{x)lpn{x) = Sm,n ■ (21) 



2m 

The Kohn-Sham energies e„ are the diagonal elements of the Lagrange multiplier A, when the 
latter is forced to be diagonal. 

The requirement that the derivative of the total energy Eq. |20] with respect to the density van- 
ishes, yields an expression for the effective potential 

J 47reo|r- r'| SnW{r) 
Both equations, together with the density constraint 



form a set of coupled equations, that determine the electron density and the total energy. This 
set of coupled equations, Eqs.[2Tl[22l and|23l is what is solved in the so-called self-consistency 
loop. Once the set of self-consistent equations has been solved, we obtain the electron density 
and we can evaluate the total energy. 
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Fig. 3: Self-consistency cycle 



In practice, one often makes the assumption that the non-interacting electrons in the effective 
potential closely resemble the true interacting electrons, and extracts a wealth of other physical 
properties from the Kohn-Sham wave functions \4'n) and the Kohn-Sham energies e.„. However, 
there is little theoretical backing for this approach and, if it fails, one should not blame density 
functional theory! 



Is there a density functional? 

The argument leading to the self-consistent equations, Eqs.[2Tl[22l and[23l relied entirely on 
the hope that exchange correlation functional can be expressed as a functional of the electron 
density. In fact, this can easily be shown, if we restrict us to ground state densities. The proof 
goes back to the seminal paper by Levy lfT5l[T6l . 

Imagine that one could construct all fermionic many-particle wave functions. For each of these 
wave functions, we can determine in a unique way the electron density 

nW(f) = N^j d^X2 ... j d^XN \^{x,X2, ...xn)\^ . (24) 

Having the electron densities, we sort the wave functions according to their density. For each 
density, I get a mug M[?7.'^^^] that holds all wave functions with that density, that is written on 
the label of the mug. 

Now we turn to each mug M[r2(^^] in sequence and determine for each the wave function with 
the lowest energy. Because the external potential energy is the same for all wave functions with 
the same density, we need to consider only the kinetic energy operator T and the operator W of 
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Fig. 4: Illustration for Levy's proof that there exists a density functional 

the electron-electron interaction, and we do not need to consider the external potential. 

= min {^\f + W\^) (25) 

\-^)eM[nW] 

F^[n^^^ is the universal density functional. It is universal in the sense that it is an intrinsic 
property of the electron gas and absolutely independent of the external potential. 
Next, we repeat the same construction as that for a universal density functional, but now we 
leave out the interaction W and consider only the kinetic energy T. 

F^\n^^^] = min (*|f I*) (26) 
The resulting functional is nothing but the kinetic energy of non-interacting electrons 

Now we can write down the total energy as functional of the density 

= F^[n^^^] + J d^rve^tir}n^^\r) (27) 

When we compare Eq. [27] with Eq. [20l we obtain an expression for the exchange correlation 
energy. 

E,,[n'-'^] = F^[n(i)(^1] - F°[n^'\f)] I d'r I d\' ^'^^'^^l^^'^O ^28) 

2 J J 47reo|r — r'\ 
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This completes the proof that the exchange correlation energy is a functional of the electron 
density. The latter was the assumption for the derivation of the set of self-consistent equations, 
Eqs.[2Tl[22l and[23]for the Kohn-Sham wave functions tpnix). 

With this, I finish the description of the theoretical basis of density-functional theory. We have 
seen that the total energy can rigorously be expressed as a functional of the density or, in prac- 
tice, as a functional of a set of one-particle wave functions, the Kohn-Sham wave functions and 
their occupations. Density functional theory per se is not an approximation and, in contrast to 
common belief, it is not a mean-field approximation. Nevertheless, we need to introduce ap- 
proximations to make density functional theory work. This is because the exchange correlation 
energy Exc[n^^^] is not completely known. These approximations will be discussed in the next 
section. 

3 Jacob's ladder of density functionals 

The development of density functionals is driven by mathematical analysis of the exact ex- 
change correlation hole [[1411711. physical insight and numerical benchmark calculations on real 
systems. The functionals evolved in steps from one functional form to another, with several 
parameterizations at each level. Perdew pictured this development by Jacob's ladder leading up 
to heaven liTTl fTl. In his analogy the different rungs of the ladder represent the different levels 
of density functionals leading to the unreachable, ultimately correct functional. 

LDA, the big surprise 

The first density functionals used in practice were based on the local-density approximation 
(LDA). The hole function for an electron at position f has been approximated by the one of a 
homogeneous electron gas with the same density as n^^\r). The exchange correlation energy 
for the homogeneous electron gas has been obtained by quantum Monte Carlo calculations ifTSl 
and analytic calculations |[T9l . The local density approximation has been generalized early to 
local spin-density approximation (LSD) |[20l . 

Truly surprising was how well the theory worked for real systems. Atomic distances could 
be determined within a few percent of the bond length and energy differences in solids were 
surprisingly good. 

This was unexpected, because the density in real materials is far from homogeneous. Gun- 
narsson and Lundquist [[2T][ explained this finding with sumrules, that are obeyed by the local 
density approximation: Firstly, the exchange correlation energy depends only on the spherical 
average of the exchange correlation hole. Of the radial hole density only the first moment con- 
tributes, while the second moment is fixed by the sum-rule that the electron density of the hole 
integrates to —1. Thus we can use 




(29) 
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where the angular brackets imply the angular average of r' — tq. This dependence on the 
hole density is rather insensitive to small changes of the hole density. Even for an atom, the 
spherically averaged exchange hole closely resembles that of the homogeneous electron gas 

Bl. 

The main deficiency of the LDA was the strong overbinding with bond energies in error by 
about one electron volt. On the one hand, this rendered LDA useless for most applications in 
chemistry. On the other hand, the problem was hardly visible in solid state physics where bonds 
are rarely broken, but rearranged so that errors cancelled. 

GGA, entering chemistry 

Being concerned about the large density variations in real materials, one tried to include the 
first terms of a Taylor expansion in the density gradients. These attempts failed miserably. The 
culprit has been a violation of the basic sum rules as pointed out by Perdew [[22| . The cure was 
a cutoff for the gradient contributions at high gradients, which lead to the class of generalized 
gradient approximations (GGA) ll23l . 

Becke ll24ll provides an intuitive description for the workings of GGA's, which I will sketch here 
in a simplified manner: Becke uses an ansatz E^c = / d^i^ A{n{r))F{x{f)) for the exchange- 
correlation energy where n{f) is the local density and x = \Vn\/n^ is a dimensionless reduced 
gradient. Do not confuse this symbol with the combined position-and-spin coordinate x. The 
function A is simply the LDA expression and F{x) is the so-called enhancement factor. The 
large-gradient limit of F{x) is obtained from a simple physical argument: 




Fig. 5: Left figure: reduced density gradient x = |Vn|/nt of a silicon atom as function of 
distance from the nucleus demonstrating that the largest reduced gradients occur in the expo- 
nential tails. Right figure: additional contribution from the gradient correction (PBE versus 
PW91 LDA) of the exchange correlation energy per electron. The figure demonstrates that the 
gradient correction stabilizes the tails of the wave function. The covalent radius of silicon is at 

1.11 A. 

Somewhat surprisingly, the reduced gradient is largest not near the nucleus but in the exponen- 
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tially decaying charge-density tails as shown in Fig. [51 For an electron that is far from an atom, 
the hole is on the atom, because a hole can only be dug where electrons are. Thus the Coulomb 

2 

interaction energy of the electron with its hole is where r is the distance of the reference 

electron from the atom. As shown in appendix El the enhancement factor can now be obtained 
by enforcing this behavior for exponentially decaying densities. 

As a result, the exchange and correlation energy per electron in the tail region of the electron 
density falls of with the inverse distance in GGA, while it has a much faster, exponential decay 
in the LDA. Thus, the the tail region is stabilized by GGA. This contribution acts like a negative 
"surface energy". 

When a bond between two atoms is broken, the surface is increased. In GGA this bond-breaking 
process is more favorable than in LDA, and, hence, the bond is weakened. Thus the GGA cures 
the overbinding error of the LDA. 

These gradient corrections greatly improved the bond energies and made density functional 
theory useful also for chemists. The most widely distributed GGA functional is the Perdew- 
Burke-Emzerhof (PBE) functional |l25l. 



Meta GGA's 

The next level of density functionals are the so-called meta GGA's ll26l l27l l28l that include not 
only the gradient of the density, but also the second derivatives of the density. These functionals 
can be reformulated so that the additional parameter is the kinetic energy density instead of the 
second density derivatives. Perdew recommends his TPSS functional [29|. 



Hybrid functionals 

Another generation of functionals are hybrid functionals [|30l ISTI . which replace some of the 
exchange energy by the exact exchange 

Er =--y ujn [ [ ^...t^umm^i)^^ oo) 

2-^^ J J 47ren|r — r'l 

m.r). IX w U| I 



where /„ and the i>n{x) are the Kohn-Sham occupations and wave functions, respectively. 
The motivation for this approach goes back to the adiabatic connection formula OH [33l |21]| 



1 r 



E^Mm = I d\ U^^[n{r)] = I dh n(f) / d\\ [ d'r' ^'[^'l (31) 
Jo J Jo ^ J 47re|r — r\ 

which expresses the exchange correlation energy as an integral of the potential energy of ex- 
change and correlation over the interaction strength. Here the interaction in the Hamiltonian is 
scaled by a factor A, leading to a A-dependent universal functional F^^[n^^'']. The interaction 
energy can be expressed by 

Jo dX 

= nn] + \l ! dh' + f dX U^f[n] (32) 

2 J J 47reo|r — r'\ Jo 
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which leads via Eq.[28]to Eq.[3TJ Using perturbation theory, the derivative of simplifies 
to the expectation (^(A)|iy|^(A)) value of the interaction, which is the potential energy of 
exchange and correlation evaluated for a many -particle wave function obtained for the specified 
given interaction strength. 

The underlying idea of the hybrid functionals is to interpolate the integrand between the end 
points. In the non-interacting limit, i.e. for A = the integrand is exactly given by the 
exact exchange energy of Eq. [30l For the full interaction, on the other hand, the LDA or GGA 
functionals are considered correctly. Thus a linear interpolation would yield 

E.. = \ {ul + f/f ) = I [e^ + f/f ) = Eg^^ + \ [e^ - i^r^) • (33) 

Depending on whether the A-dependence is a straight line or whether it is convex, the weight 
factor may be equal or smaller than |. Perdew [|34ll has given arguments that a factor | would 
actually be better than a factor |. 

Hybrid functionals perform substantially better than GGA functionals regarding binding ener- 
gies, band gaps and reaction energies. However, they are flawed for the description of solids. 
The reason is that the exact exchange hole in a solid is very extended. These long-range tails 
are screened away quickly when the interaction is turned on, because they are cancelled by the 
correlation. Effectively, we should use a smaller mixing factor for the long range part of the 
exchange hole. This can be taken into account, by cutting off the long-range part of the interac- 
tion for the calculation of the Hartree-Fock exchange [[35l . This approach improves the results 
for band gaps while reducing the computational effort ll36l . 

The effective cancellation of the long-ranged contribution of exchange with a similar contri- 
bution from correlation, which is also considered properly already in the LDA, is one of the 
explanation for the superiority of the LDA over the Hartree-Fock approximation. 
The most widely used hybrid functional is the B3LYP functional ll37l . which is, however, ob- 
tained from a parameter fit to a database of simple molecules. The functional PBEO [[38l[39l is 
bom out of the famous PBE GGA functional and is a widely distributed parameter-free func- 
tional. 

LDA+U and local hybrid functionals 

Starting from a completely different context, Anisimov et. al. [|40l introduced the so-called 
LDA-i-U method, which, as described below, has some similarities to the hybrid functionals 
above. 

The main goal was to arrive at a proper description of transition metal oxides, which tend to 
be Mott insulators, while GGA calculations predict them often to be metals. The remedy was 
to add a correlation tenr§ [|4TI borrowed from the Hubbard model and to correct the resulting 



^The expression given here looks unusually simple. This is due to the notation of spin orbitals, which takes 
care of the spin indices. 
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double counting of the interactions by Edc- 



~'~ 2 ^ ^ f^Q,/3,7,5 f P7,aP<5,/3 — P5,qP7,/3 j — -E'dc (34) 
R a,l3,'y,5eCR ^ ^ 

t^Q,/3,7,(5 = a X a X — ^ (35) 

J J 47reo|r — r'\ 

|^n)/n(^nk/3) (36) 

where \xa) are atomic tight-binding orbitals and IvTq,) are their projector functions]^ The addi- 
tional energy is a Hartree-Fock exchange energy, that only considers the exchange for specified 
sets of local orbitals. The exchange term does only consider a subset of orbitals Cr for each 
atom R and it ignores the contribution involving orbitals centered on different atoms. 
Novak et al. Il42l made the connection to the hybrid functionals explicit and restricted the 
exact exchange contribution of a hybrid functional to only a shell of orbitals. While in the 
LDA-i-U method the bare Coulomb matrix elements are reduced by a screening factor, in the 
hybrid functionals it is the mixing factor that effectively plays the same role. Both, LDA-i-U 
and the local hybrid method, have in common that they radically remove the contribution of 
off-site matrix elements of the interaction. Tran et al. H3l applied this method to transition 
metal oxides and found results that are similar to those of the full implementation of hybrid 
functionals. 



Van der Waals interactions 

One of the major difficulties for density functionals is the description of van der Waals forces, 
because it is due to the quantum mechanical synchronization of charge fluctuations on distinct 
molecules. I refer the reader to the work made in the group of Lundqvist p4ll45ll46l . 



4 Benchmarks, successes and failures 

The development of density functionals has profited enormously on careful benchmark studies. 
The precondition is a data set of test cases for which reliable and accurate experimental data 
exist. The most famous data sets are the Gl and G2 databases ll47l l48l |49l l50l that have been 
set up to benchmark quantum-chemistry codes. Becke lISTl [52l |3T1 l53l [54l set a trend by using 
these large sets of test cases for systematic studies of density functionals. In order to separate 
out the accuracy of the density functionals, it is vital to perform these calculations on extremely 
accurate numerical methods. Becke used basis set free calculations, that were limited to small 
molecules while being extremely accurate. Paier et. al. Il55l l56l 1571 [36l have later performed 

^Projector functions obey the biorthogonality condition (xqItt^) = 5a,i3- Within the sub-Hilbert space of the 
tight-binding orbitals, i.e. for wave functions of the form \'4>) = |Xq)cq, the projector functions decompose 
the wave function into tight binding orbitals, i.e. j^/') ~ (""q IV^)- A similar projection is used extensively 

in the projector augmented-wave method described later 
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careful comparisons of two methods, Gaussian and the projector augmented-wave method, to 
single out the error of the electronic structure method. 

Overall, the available density functionals predict molecular structures very well. Bond distances 
agree with the experiment often within one percent. Bond angles come out within a few degrees. 
The quality of total energies depends strongly on the level of functionals used. On the LDA 
level bonds are overestimated in the 1 eV range, on the GGA level these errors are reduced to 
a about 0.3 eV, and hybrid functionals reduce the error by another factor of 2. The ultimate 
goal is to reach chemical accuracy, which is about 0.05 eV. Such an accuracy allows to predict 
reaction rates at room temperature within a factor of 10. 

Band gaps are predicted too small with LDA and GGA. The so-called band gap problem has 
been one of the major issues during the development of density functionals. Hybrid functionals 
clearly improve the situation. A problem is the description of materials with strong electron 
correlations. For LDA and GGA many insulating transition metal oxides are described as met- 
als. This changes again for the hybrid functionals, which turns them into antiferromagnetic 
insulators, which is a dramatic improvement. 



5 Electronic structure methods 



In this second part of my lecture notes, I will address the problem how to solve the Kohn-Sham 
equations and how to obtain the total energy and other observables. It is convenient to use a 
slightly different notation: Instead of treating the nuclei via an external potential, we combine 
all electrostatic interactions into a single double integral. 
This brings the total energy into the form 



E 



{^„(f)},{i?H} 
1 

+5 



p ■ 



n{r') + Z{r' 



47ren|r — r' 



(37) 



where Z{r) = — ZR6{f— Rr) is the nuclear charge density expressed in electron charges. 
Zji is the atomic number of a nucleus at position Rji. 

The electronic ground state is determined by minimizing the total energy functional of 
Eq. |37]at a fixed ionic geometry. The one-particle wave functions have to be orthogonal. This 
constraint is implemented with the method of Lagrange multipliers. We obtain the ground-state 
wave functions from the extremum condition for 



Y 



E 



Ar 



(38) 



with respect to the wavefunctions and the Lagrange multipliers 
for the wavefunctions has the form 



The extremum condition 



H\lpn)fn = ^ \^m)Am,r 



(39) 
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where H = i^p"^ + ?)cff is the effective one-particle Hamilton operator. 

The corresponding effective potential depends itself on the electron density via 

p e'^(n{r') + Z{r')] 

^effi^ = / d'r' ^, ^ + /i.c(rO , (40) 

J 47reo|r — r'l 

where jixcir) = ^^I'^pp'^ is the functional derivative of the exchange and correlation functional. 
After a unitary transformation that diagonalizes the matrix of Lagrange multipliers A, we obtain 
the Kohn-Sham equations 

H\ijn) = men . (41) 
The one-particle energies e„ are the eigenvalues of the matrix with the elements Is.n,m{fn + 

fm)/{2fnfm) ffSl. 

The one-electron Schrodinger equations, namely the Kohn-Sham equations given in Eq. [211 
still pose substantial numerical difficulties: (1) in the atomic region near the nucleus, the kinetic 
energy of the electrons is large, resulting in rapid oscillations of the wavefunction that require 
fine grids for an accurate numerical representation. On the other hand, the large kinetic energy 
makes the Schrodinger equation stiff, so that a change of the chemical environment has little 
effect on the shape of the wavefunction. Therefore, the wavefunction in the atomic region can 
be represented well already by a small basis set. (2) In the bonding region between the atoms 
the situation is opposite. The kinetic energy is small and the wavefunction is smooth. However, 
the wavefunction is flexible and responds strongly to the environment. This requires large and 
nearly complete basis sets. 

Combining these different requirements is non-trivial and various strategies have been devel- 
oped. 

• The atomic point of view has been most appealing to quantum chemists. Basis functions 
are chosen that resemble atomic orbitals. This choice exploits that the wavefunction in 
the atomic region can be described by a few basis functions, while the chemical bond is 
described by the overlapping tails of these atomic orbitals. Most techniques in this class 
are a compromise of, on the one hand, a well adapted basis set, where the basis functions 
are difficult to handle, and, on the other hand, numerically convenient basis functions 
such as Gaussians, where the inadequacies are compensated by larger basis sets. 

• Pseudopotentials regard an atom as a perturbation of the free electron gas. The most nat- 
ural basis functions for the free electron gas are plane waves. Plane-wave basis sets are in 
principle complete and suitable for sufficiently smooth wavefunctions. The disadvantage 
of the comparably large basis sets required is offset by their extreme numerical simplicity. 
Finite plane-wave expansions are, however, absolutely inadequate to describe the strong 
oscillations of the wavefunctions near the nucleus. In the pseudopotential approach the 
Pauli repulsion by the core electrons is therefore described by an effective potential that 
expels the valence electrons from the core region. The resulting wavefunctions are smooth 
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and can be represented well by plane waves. The price to pay is that all information on 
the charge density and wavefunctions near the nucleus is lost. 

• Augmented-wave methods compose their basis functions from atom-like wavefunctions 
in the atomic regions and a set of functions, called envelope functions, appropriate for the 
bonding in between. Space is divided accordingly into atom-centered spheres, defining 
the atomic regions, and an interstitial region in between. The partial solutions of the 
different regions, are matched with value and derivative at the interface between atomic 
and interstitial regions. 

The projector augmented-wave method is an extension of augmented wave methods and the 
pseudopotential approach, which combines their traditions into a unified electronic structure 
method. 

After describing the underlying ideas of the various approaches, let us briefly review the history 
of augmented wave methods and the pseudopotential approach. We do not discuss the atomic- 
orbital based methods, because our focus is the PAW method and its ancestors. 

6 Augmented wave methods 

The augmented wave methods have been introduced in 1937 by Slater ll59l . His method was 
called augmented plane- wave (APW) method. Later Korringa ll60l . Kohn and Rostokker ||6TI 
modified the idea, which lead to the so-called KKR method. The basic idea behind the aug- 
mented wave methods has been to consider the electronic structure as a scattered-electron prob- 
lem: Consider an electron beam, represented by a plane wave, traveling through a solid. It 
undergoes multiple scattering at the atoms. If, for some energy, the outgoing scattered waves 
interfere destructively, so that the electrons can not escape, a bound state has been determined. 
This approach can be translated into a basis-set method with energy- and potential-dependent 
basis functions. In order to make the scattered wave problem tractable, a model potential had 
to be chosen: The so-called muffin-tin potential approximates the true potential by a potential, 
that is spherically symmetric in the atomic regions, and constant in between. 
Augmented-wave methods reached adulthood in the 1970s: O. K. Andersen |[62l showed that 
the energy dependent basis set of Slater's APW method can be mapped onto one with energy in- 
dependent basis functions by linearizing the partial waves for the atomic regions with respect to 
their energy. In the original APW approach, one had to determine the zeros of the determinant 
of an energy dependent matrix, a nearly intractable numerical problem for complex systems. 
With the new energy independent basis functions, however, the problem is reduced to the much 
simpler generalized eigenvalue problem, which can be solved using efficient numerical tech- 
niques. Furthermore, the introduction of well defined basis sets paved the way for full-potential 
calculations [633. In that case, the muffin-tin approximation is used solely to define the basis 
set Ixi), while the matrix elements {xi\H\xj) of the Hamiltonian are evaluated with the full 
potential. 
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In the augmented wave methods one constructs the basis set for the atomic region by solving 
the radial Schrodinger equation for the spherically averaged effective potential 



2m, 



+ V. 



eff\ 







as function of the energy. Note that a partial wave 0£,m(e, r) is an angular-momentum eigenstate 
and can be expressed as a product of a radial function and a spherical harmonic. The energy- 
dependent partial wave is expanded in a Taylor expansion about some reference energy e^^^ 



where 4)u,t,m{r) = 4>i,m{^u,t-, r). The energy derivative of the partial wave (p^ir) 
obtained from the energy derivative of the Schrodinger equation 
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Next, one starts from a regular basis set, such as plane waves, Gaussians or Hankel functions. 
These basis functions are called envelope functions Within the atomic region they are 
replaced by the partial waves and their energy derivatives, such that the resulting wavefunction 
Xi{r) is continuous and differentiable. The augmented envelope function has the form 



R 



R,l,m 



'Pu,R,e,m{f^O,R,i,m,i + 4>u,R,£,m ,m,i 



(42) 



OR{f) is a step function that is unity within the augmentation sphere centered at Rr and zero 
elsewhere. The augmentation sphere is atom-centered and has a radius about equal to the cova- 
lent radius. This radius is called the muffin-tin radius, if the spheres of neighboring atoms touch. 
These basis functions describe only the valence states; the core states are localized within the 
augmentation sphere and are obtained directly by a radial integration of the Schrodinger equa- 
tion within the augmentation sphere. 

The coefficients aji^e,m,i and bji^i^rn,i are obtained for each as follows: The envelope function 
is decomposed around each atomic site into spherical harmonics multiplied by radial functions 



Xi{r) = E?^i?,f,m,i(|r - RnDYe^rnir- Rr) . 



(43) 



£,m 



Analytical expansions for plane waves, Hankel functions or Gaussians exist. The radial parts 
of the partial waves (f)u,R,e,m and (j)u,R,e,m are matched with value and derivative to MK/,m,j(|^), 
which yields the expansion coefficients aR^i^m,i and hR/^m,i- 

If the envelope functions are plane waves, the resulting method is called the linear augmented 
plane-wave (LAPW) method. If the envelope functions are Hankel functions, the method is 
called linear muffin-tin orbital (LMTO) method. 

A good review of the LAPW method [[62l has been given by Singh ll64ll . Let us now briefly men- 
tion the major developments of the LAPW method: Soler ll65l introduced the idea of additive 
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augmentation: While augmented plane waves are discontinuous at the surface of the augmen- 
tation sphere if the expansion in spherical harmonics in Eq.|42]is truncated, Soler replaced the 
second term in Eq. |42l by an expansion of the plane wave with the same angular momentum 
truncation as in the third term. This dramatically improved the convergence of the angular mo- 
mentum expansion. Singh [[66l introduced so-called local orbitals, which are non-zero only 
within a muffin-tin sphere, where they are superpositions of and functions from different 
expansion energies. Local orbitals substantially increase the energy transferability. Sjostedt 
ll67l relaxed the condition that the basis functions are differentiable at the sphere radius. In 
addition she introduced local orbitals, which are confined inside the sphere, and that also have 
a kink at the sphere boundary. Due to the large energy cost of kinks, they will cancel, once the 
total energy is minimized. The increased variational degree of freedom in the basis leads to a 
dramatically improved plane- wave convergence 

The second variant of the linear methods is the LMTO method [[62l . A good introduction into 
the LMTO method is the book by Skriver ll69ll . The LMTO method uses Hankel functions as 
envelope functions. The atomic spheres approximation (ASA) provides a particularly simple 
and efficient approach to the electronic structure of very large systems. In the ASA the aug- 
mentation spheres are blown up so that the sum of their volumes is equal to the total volume. 
Then, the first two terms in Eq.|42]are ignored. The main deficiency of the LMTO- ASA method 
is the limitation to structures that can be converted into a closed packed arrangement of atomic 
and empty spheres. Furthermore, energy differences due to structural distortions are often qual- 
itatively incorrect. Full potential versions of the LMTO method, that avoid these deficiencies 
of the ASA have been developed. The construction of tight binding orbitals as superposition of 
muffin-tin orbitals [iTOl showed the underlying principles of the empirical tight-binding method 
and prepared the ground for electronic structure methods that scale linearly instead of with the 
third power of the number of atoms. The third generation LMTO ITtTI allows to construct true 
minimal basis sets, which require only one orbital per electron pair for insulators. In addition, 
they can be made arbitrarily accurate in the valence band region, so that a matrix diagonalization 
becomes unnecessary. The first steps towards a full-potential implementation, that promises a 
good accuracy, while maintaining the simplicity of the LMTO-ASA method are currently under 
way. Through the minimal basis-set construction the LMTO method offers unrivaled tools for 
the analysis of the electronic structure and has been extensively used in hybrid methods com- 
bining density functional theory with model Hamiltonians for materials with strong electron 
correlations ir72l . 

7 Pseudopotentials 

Pseudopotentials have been introduced to (1) avoid describing the core electrons explicitly and 
(2) to avoid the rapid oscillations of the wavefunction near the nucleus, which normally require 
either complicated or large basis sets. 

The pseudopotential approach traces back to 1940 when C. Herring invented the orthogonalized 
plane-wave method fl3\. Later, Phillips ll74ll and Antoncik [[751 replaced the orthogonality 
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condition by an effective potential, which mimics the Pauli repulsion by the core electrons 
and thus compensates the electrostatic attraction by the nucleus. In practice, the potential was 
modified, for example, by cutting off the singular potential of the nucleus at a certain value. This 
was done with a few parameters that have been adjusted to reproduce the measured electronic 
band structure of the corresponding solid. 

Hamann, Schliiter and Chiang [1761 showed in 1979 how pseudopotentials can be constructed 
in such a way, that their scattering properties are identical to that of an atom to first order in 
energy. These first-principles pseudopotentials relieved the calculations from the restrictions of 
empirical parameters. Highly accurate calculations have become possible especially for semi- 
conductors and simple metals. An alternative approach towards first-principles pseudopotentials 
by Zunger and CohenfTT] even preceded the one mentioned above. 



The idea behind the pseudopotential construction 

In order to construct a first-principles pseudopotential, one starts out with an all-electron density- 
functional calculation for a spherical atom. Such calculations can be performed efficiently on 
radial grids. They yield the atomic potential and wavefunctions </)^,m(r). Due to the spherical 
symmetry, the radial parts of the wavefunctions for different magnetic quantum numbers m are 
identical. 

For the valence wavefunctions one constructs pseudo wavefunctions |0^,m) • There are numerous 
ways iTTSl 1791 [80l [STIl to construct those pseudo wavefunctions: Pseudo wave functions are 
identical to the true wave functions outside the augmentation region, which is called core region 
in the context of the pseudopotential approach. Inside the augmentation region the pseudo 
wavefunction should be node-less and have the same norm as the true wavefunctions, that is 
(0£,m|0£,m) = {<Pi,m\<pLm) (comparc Figure©. 

From the pseudo wavefunction, a potential ui (r ) can be reconstructed by inverting the respective 
Schrodinger equation, i.e. 

1 ^ ~ 

(t>i,m{r) = ^ ui{r) = e + ^ — ■ - — V>£,m(f) . 

This potential u^if) (compare Figure [6l), which is also spherically symmetric, differs from one 
main angular momentum l to the other. Note, that this inversion of the Schrodinger equation 
works only if the wave functions are nodeless. 
Next we define an effective pseudo Hamiltonian 

Hi = -T^V^ + (f) + / d'r'^ ^, ^ + yU.c([n(r-0], r) , (44) 

2me J 47reo|f-r'| 

where fixc{r) = SExc[n]/Sn{f) is the functional derivative of the exchange and correlation 
energy with respect to the electron density. Then, we determine the pseudopotentials v^'^ such 
that the pseudo Hamiltonian produces the pseudo wavefunctions, that is 

r e^(fi{f') + Z{f')) 

vf{r) = M,(f) - / ^ ,^ ^, ^ - /i,.e([n(F)],F) . (45) 

J Aireolr — r'\ 
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Fig. 6: Illustration of the pseudopotential concept at the example of the 3s wavefunction of Si. 
The solid line shows the radial part of the pseudo wavefunction (pi^m- The dashed line corre- 
sponds to the all-electron wavefunction (pi^m^ which exhibits strong oscillations at small radii. 
The angular momentum dependent pseudopotential ui (dash-dotted line) deviates from the all- 
electron potential Veff (dotted line) inside the augmentation region. The data are generated by 
the fhi98PP code W\ . 



This process is called "unscreening". 

Z{f) mimics the charge density of the nucleus and the core electrons. It is usually an atom- 
centered, spherical Gaussian that is normalized to the charge of nucleus and core of that atom. 
In the pseudopotential approach, ZR{r) does not change with the potential. The pseudo density 

= Ylin fni^ni^)i'n{r) is constructcd from the pseudo wavefunctions. 
In this way, we obtain a different potential for each angular momentum channel. In order to 
apply these potentials to a given wavefunction, the wavefunction must first be decomposed 
into angular momenta. Then each component is applied to the pseudopotential v^'^ for the 
corresponding angular momentum. 

The pseudopotential defined in this way can be expressed in a semi-local form 

yP^(^f,r') = v{f)5{r- r^) + 

The local potential v{f) only acts on those angular momentum components that are not already 
considered explicitely in the non-local, angular-momentum dependend pseudopotentials f J*. 
Typically it is chosen to cancel the most expensive nonlocal terms, the one corresponding to the 
highest physically relevant angular momentum. 

The pseudopotential v^^{f, r') is non-local as its depends on two position arguments, r and f" . 
The expectation values are evaluated as a double integral 

{i)\vps\i') = I d^r I d^r' ij*{r)vP'{fy)ip{f') (47) 



Y(>,mir) [ff (r) - v{r}] 
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The semi-local form of the pseudopotential given in Eq. |46] is computationally expensive. 
Therefore, in practice one uses a separable form of the pseudopotential [[83l [84l [85l 

~ . (48) 



Thus, the projection onto spherical harmonics used in the semi-local form of Eq.|46]is replaced 
by a projection onto angular momentum dependent functions vP^\(j)i). 

The indices i and j are composite indices containing the atomic-site index R, the angular mo- 
mentum quantum numbers m and an additional index a. The index a distinguishes partial 
waves with otherwise identical indices i?, £, m when more than one partial wave per site and 
angular momentum is allowed. The partial waves may be constructed as eigenstates to the 
pseudopotential v^/ for a set of energies. 

One can show that the identity of Eq.|48] holds by applying a wavefunction = J2i to 
both sides. If the set of pseudo partial waves in Eq. |48]is complete, the identity is exact. 
The advantage of the separable form is that {(f)\v'^'^ is treated as one function, so that expectation 
values are reduced to combinations of simple scalar products {(pilv^^lip). 
The total energy of the pseudopotential method can be written in the form 



n I ^ps 



2m^ 

n 



1 r r e^fhir) + Z{r))(n{r') + Z{f 



+- / d'r / d'r'^ + E.,[n(f)] . (49) 

The constant Egeij is adjusted such that the total energy of the atom is the same for an all- 
electron calculation and the pseudopotential calculation. 

For the atom, from which it has been constructed, this construction guarantees that the pseu- 
dopotential method produces the correct one-particle energies for the valence states and that the 
wave functions have the desired shape. 

While pseudopotentials have proven to be accurate for a large variety of systems, there is no 
strict guarantee that they produce the same results as an all-electron calculation, if they are used 
in a molecule or solid. The error sources can be divided into two classes: 

• Energy transferability problems: Even for the potential of the reference atom, the scatter- 
ing properties are accurate only in given energy window. 

• Charge transferability problems: In a molecule or crystal, the potential differs from that 
of the isolated atom. The pseudopotential, however, is strictly valid only for the isolated 
atom. 

The plane-wave basis set for the pseudo wavefunctions is defined by the shortest wave length 
A = 2n/\G\, where G is the wave vector, via the so-called plane- wave cutoff Epw = ^ 
with Gmax = max{|G|}. It is often specified in Rydberg (1 Ry=| H^13.6 eV). The plane- 
wave cutoff is the highest kinetic energy of all basis functions. The basis-set convergence can 
systematically be controlled by increasing the plane- wave cutoff. 
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The charge transferability is substantially improved by including a nonlinear core correction 
lf86l into the exchange-correlation term of Eq. |49l Hamann ifSTl showed, how to construct 
pseudopotentials also from unbound wavefunctions. Vanderbilt ll88l[89l generalized the pseu- 
dopotential method to non-normconserving pseudopotentials, so-called ultra-soft pseudopoten- 
tials, which dramatically improves the basis-set convergence. The formulation of ultra-soft 
pseudopotentials has already many similarities with the projector augmented-wave method. 
Truncated separable pseudopotentials suffer sometimes from so-called ghost states. These are 
unphysical core-like states, which render the pseudopotential useless. These problems have 
been discussed by Gonze |[90l . Quantities such as hyperfine parameters that depend on the full 
wavefunctions near the nucleus, can be extracted approximately [|9T]| . A good review about 
pseudopotential methodology has been written by Payne ll92l and Singh [|64ll . 

In 1985 R. Car and M. Parrinello published the ab-initio molecular dynamics method [|93l . 
Simulations of the atomic motion have become possible on the basis of state-of-the-art elec- 
tronic structure methods. Besides making dynamical phenomena and finite temperature effects 
accessible to electronic structure calculations, the ab-initio molecular dynamics method also 
introduced a radically new way of thinking into electronic structure methods. Diagonalization 
of a Hamilton matrix has been replaced by classical equations of motion for the wavefunction 
coefficients. If one applies friction, the system is quenched to the ground state. Without fric- 
tion truly dynamical simulations of the atomic structure are performed. By using thermostats 
ll94l l95l l96l l97l . simulations at constant temperature can be performed. The Car-Parrinello 
method treats electronic wavefunctions and atomic positions on an equal footing. 



8 Projector augmented-wave method 

The Car-Parrinello method had been implemented first for the pseudopotential approach. There 
seemed to be insurmountable barriers against combining the new technique with augmented 
wave methods. The main problem was related to the potential dependent basis set used in 
augmented wave methods: the Car-Parrinello method requires a well defined and unique to- 
tal energy functional of atomic positions and basis set coefficients. Furthermore the analytic 
evaluation of the first partial derivatives of the total energy with respect to wave functions, 
= H\'ijjn)fn, and atomic positions, the forces Fj = —VjE, must be possible. Therefore, 
it was one of the main goals of the PAW method to introduce energy and potential independent 
basis sets, which were as accurate as the previously used augmented basis sets. Other require- 
ments have been: (1) The method should at least match the efficiency of the pseudopotential 
approach for Car-Parrinello simulations. (2) It should become an exact theory when converged 
and (3) its convergence should be easily controlled. We believe that these criteria have been 
met, which explains why the PAW method becomes increasingly wide spread today. 
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Transformation theory 

At the root of the PAW method lies a transformation, that maps the true wavefunctions with their 
complete nodal structure onto auxiliary wavefunctions, that are numerically convenient. We aim 
for smooth auxiliary wavefunctions, which have a rapidly convergent plane-wave expansion. 
With such a transformation we can expand the auxiliary wave functions into a convenient basis 
set such as plane waves, and evaluate all physical properties after reconstructing the related 
physical (true) wavefunctions. 

Let us denote the physical one-particle wavefunctions as |^„) and the auxiliary wavefunctions 
as IV^n)- Note that the tilde refers to the representation of smooth auxiliary wavefunctions and 
n is the label for a one-particle state and contains a band index, a /c-point and a spin index. The 
transformation from the auxiliary to the physical wave functions is denoted by T, i.e. 

IV'n) = r|^„) . (50) 



Now we express the constrained density functional F of Eq. [38]in terms of our auxiliary wave- 
functions 



Am,n • (51) 



The variational principle with respect to the auxiliary wavefunctions yields 

pHf\tPn)=Pf\i^n)en. (52) 

Again, we obtain a Schrodinger-like equation (see derivation of Eq.l4T)). but now the Hamilton 
operator has a different form, H = T^HT, an overlap operator O = T'^T occurs, and the 
resulting auxiliary wavefunctions are smooth. 

When we evaluate physical quantities, we need to evaluate expectation values of an operator A, 
which can be expressed in terms of either the true or the auxiliary wavefunctions, i.e. 

(i) = 5^/n(^n|i|V'n) = $^/n(^An|rUr|^/;„). (53) 
n n 

In the representation of auxiliary wavefunctions we need to use transformed operators A = 
AT. As it is, this equation only holds for the valence electrons. The core electrons are 
treated differently as will be shown below. 

The transformation takes us conceptionally from the world of pseudopotentials to that of aug- 
mented wave methods, which deal with the full wavefunctions. We will see that our auxiliary 
wavefunctions, which are simply the plane-wave parts of the full wavefunctions, translate into 
the wavefunctions of the pseudopotential approach. In the PAW method the auxiliary wavefunc- 
tions are used to construct the true wavefunctions and the total energy functional is evaluated 
from the latter. Thus it provides the missing link between augmented wave methods and the 
pseudopotential method, which can be derived as a well-defined approximation of the PAW 
method. 
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In the original paper [58 j, the auxiliary wavefunctions have been termed pseudo wavefunctions 
and the true wavefunctions have been termed all-electron wavefunctions, in order to make the 
connection more evident. We avoid this notation here, because it resulted in confusion in cases, 
where the correspondence is not clear-cut. 

Transformation operator 

So far we have described how we can determine the auxiliary wave functions of the ground 
state and how to obtain physical information from them. What is missing is a definition of the 
transformation operator T. 

The operator T has to modify the smooth auxiliary wave function in each atomic region, so that 
the resulting wavefunction has the correct nodal structure. Therefore, it makes sense to write 
the transformation as identity plus a sum of atomic contributions Sr 

r=i + 5^5R. (54) 

R 

For every atom, Sr adds the difference between the true and the auxiliary wavefunction. 
The local terms Sr are defined in terms of solutions of the Schrodinger equation for the 
isolated atoms. This set of partial waves |</)j) will serve as a basis set so that, near the nucleus, 
all relevant valence wavefunctions can be expressed as superposition of the partial waves with 
yet unknown coefficients as 

V'(^) = for |f - Rr\ < r^^R . (55) 

ieR 

With z G -R we indicate those partial waves that belong to site R. 

Since the core wavefunctions do not spread out into the neighboring atoms, we will treat them 
differently. Currently we use the frozen-core approximation, which imports the density and 
the energy of the core electrons from the corresponding isolated atoms. The transformation T 
shall produce only wavefunctions orthogonal to the core electrons, while the core electrons are 
treated separately. Therefore, the set of atomic partial waves |</)j) includes only valence states 
that are orthogonal to the core wavefunctions of the atom. 

For each of the partial waves we choose an auxiliary partial wave The identity 

10,) = {I + SrM) for teR 
Sr^ = 10.) -10.) (56) 

defines the local contribution Sr to the transformation operator. Since 1 + iS/j shall change 
the wavefunction only locally, we require that the partial waves |0i) and their auxiliary counter 
parts |0j) are pairwise identical beyond a certain radius r^/j. 

Mr) = Mr) for i e R and \f-RR\ > r^^R (57) 

Note that the partial waves are not necessarily bound states and are therefore not normalizable, 
unless we truncate them beyond a certain radius r^i?. The PAW method is formulated such that 
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the final results do not depend on the location where the partial waves are truncated, as long 
as this is not done too close to the nucleus and identical for auxiliary and all-electron partial 
waves. 

In order to be able to apply the transformation operator to an arbitrary auxiliary wavefunction, 
we need to be able to expand the auxiliary wavefunction locally into the auxiliary partial waves 

i'ir) = ^0i(f)Q = ^0i(r)(pi|V^) for |f - Rr\ < r^^R , (58) 

which defines the projector functions \pi). The projector functions probe the local character of 
the auxiliary wave function in the atomic region. Examples of projector functions are shown in 
Figure |71 From Eq.[58]we can derive J2ieR. ~ which is valid within Tc r. It can be 

shown by insertion, that the identity Eq.[5H] holds for any auxiliary wavefunction [tp) that can 
be expanded locally into auxiliary partial waves if 

{Pi\(pj) = kj for hj^R- (59) 

Note that neither the projector functions nor the partial waves need to be orthogonal among 
themselves. The projector functions are fully determined with the above conditions and a clo- 
sure relation, which is related to the unscreening of the pseudopotentials (see Eq. 90 in ll58l ). 



Fig. 7: Projector functions of the chlorine atom. Top: two s-type projector functions, middle: 
p-type, bottom: d-type. 

By combining Eq.[56]and Eq.[58l we can apply Sr to any auxiliary wavefunction. 

SrI^p) = $^5i?|0.)(p#) = 5^(|0,)-|0i))(p,K'). (60) 

i&R i&R 

Hence, the transformation operator is 

f=l + J2{\<P^)-\4>^)){P^\, (61) 
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where the sum runs over all partial waves of all atoms. The true wave function can be expressed 



as 



i R 

with 

ieR 

li'R) = (64) 

ieR 

In Fig.[8]the decomposition of Eq.[62]is shown for the example of the bonding p-a state of the 
CI2 molecule. 




Fig. 8: Bonding p-a orbital of the CI2 molecule and its decomposition of the wavefunction into 
auxiliary wavefunction and the two one-center expansions. Top-left: True and auxiliary wave 
function; top-right: auxiliary wavefunction and its partial wave expansion; bottom-left: the two 
partial wave expansions; bottom-right: true wavefunction and its partial wave expansion. 

To understand the expression Eq.[62]for the true wave function, let us concentrate on different 
regions in space. (1) Far from the atoms, the partial waves are, according to Eq. [57l pairwise 
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identical so that the auxiliary wavefunction is identical to the true wavefunction, that is ip{r) = 
ip(r). (2) Close to an atom R, however, the auxiliary wavefunction is, according to Eg. [58l 
identical to its one-center expansion, that is ijj{r) = V'fj(^)- Hence the true wavefunction tp^f) 
is identical to which is built up from partial waves that contain the proper nodal structure. 

In practice, the partial wave expansions are truncated. Therefore, the identity of Eq. [58] does 
not hold strictly. As a result, the plane waves also contribute to the true wavefunction inside 
the atomic region. This has the advantage that the missing terms in a truncated partial wave 
expansion are partly accounted for by plane waves. This explains the rapid convergence of the 
partial wave expansions. This idea is related to the additive augmentation of the LAPW method 
of Soler Il65l. 

Frequently, the question comes up, whether the transformation Eq.[6l]of the auxiliary wavefunc- 
tions indeed provides the true wavefunction. The transformation should be considered merely 
as a change of representation analogous to a coordinate transform. If the total energy functional 
is transformed consistently, its minimum will yield auxiliary wavefunctions that produce the 
correct wave functions lip). 

Expectation values 

Expectation values can be obtained either from the reconstructed true wavefunctions or directly 
from the auxiliary wave functions 

(A) = j2ui;^\Am + 

n n=l 

= 5^/„(^„|rUr|^„) + 5^R|i|C), (65) 

n n=l 

where /„ are the occupations of the valence states and Nc is the number of core states. The first 
sum runs over the valence states, and second over the core states \(f)n)- 

Now we can decompose the matrix element for a wavefunction ^ into its individual contribu- 
tions according to Eq. [62l 

R R' 
R 

' V ' 

part 1 

R 

V ' 

part 2 

Rf^R' 

' V ' 

part 3 
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Only the first part of Eq. [66l is evaluated explicitly, while the second and third parts of Eq. [66] 
are neglected, because they vanish for sufficiently local operators as long as the partial wave 
expansion is converged: The function ijj}^ — vanishes per construction beyond its augmen- 
tation region, because the partial waves are pairwise identical beyond that region. The function 
^ ~ vanishes inside its augmentation region, if the partial wave expansion is sufficiently 
converged. In no region of space both functions tp}^ — tp\ and ^ — are simultaneously 
nonzero. Similarly the functions ^jj]^ — ip}^ from different sites are never non-zero in the same 
region in space. Hence, the second and third parts of Eq. [66] vanish for operators such as the 
kinetic energy ^^V^ and the real space projection operator |r) (r|, which produces the electron 
density. For truly nonlocal operators the parts 2 and 3 of Eq. [66] would have to be considered 
explicitly. 

The expression, Eq. [65l for the expectation value can therefore be written with the help of Eq. [66] 
as 

n n=l 

Nc ^ 

n n=l 

Nc,R 

if. ijeif nGif 

Nc,R 

R ijeR neR 
where D is the one-center density matrix defined as 

^iJ = ^fn{i'n\Pj){Pt\i'n) = ^{Pili'n) fn{i>n\Pj) ■ (68) 
n n 

The auxiliary core states, allow to incorporate the tails of the core wavefunction into the 
plane-wave part, and therefore assure, that the integrations of partial wave contributions cancel 
each other strictly beyond Tc. They are identical to the true core states in the tails, but are a 
smooth continuation inside the atomic sphere. It is not required that the auxiliary wave functions 
are normalized. 

Following this scheme, the electron density is given by 

n{r) = n{r) + ''^^(^n]^{r) — h]^{r)^ (69) 

R 

n 

i,jeR 

nUr) = 5^ Aj0*(r)0,(r-')+n,,«(f), (70) 

ijGif 
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where ric^ji is the core density of the corresponding atom and fic,K is the auxiliary core density 
that is identical to Uc^ji outside the atomic region, but smooth inside. 

Before we continue, let us discuss a special point: The matrix elements of a general operator 
with the auxiliary wavefunctions may be slowly converging with the plane- wave expansion, be- 
cause the operator A may not be well behaved. An example for such an operator is the singular 
electrostatic potential of a nucleus. This problem can be alleviated by adding an "intelligent 
zero": If an operator B is purely localized within an atomic region, we can use the identity 
between the auxiliary wavefunction and its own partial wave expansion 

= ii^nimn) - ii'lmn) ■ (71) 



Now we choose an operator 13 so that it cancels the problematic behavior of the operator A, 
but is localized in a single atomic region. By adding B to the plane-wave part and the matrix 
elements with its one-center expansions, the plane- wave convergence can be improved without 
affecting the converged result. A term of this type, namely v will be introduced in the next 
section to cancel the Coulomb singularity of the potential at the nucleus. 



Total energy 

Like wavefunctions and expectation values, also the total energy can be divided into three parts. 



E 



{\^n)}ARR}] = e + J2{e'r-e},) (72) 

R 



The plane-wave part E involves only smooth functions and is evaluated on equi-spaced grids 
in real and reciprocal space. This part is computationally most demanding, and is similar to the 
expressions in the pseudopotential approach. 



4:7ieo\f— P\ 

+ I £r v{r)h{r) + E^c[n\ (73) 



Z{r) is an angular- momentum dependent core-like density that will be described in detail below. 
The remaining parts can be evaluated on radial grids in a spherical-harmonics expansion. The 
nodal structure of the wavefunctions can be properly described on a logarithmic radial grid that 
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becomes very fine near the nucleus, 



ijGR n&R 

+ ^ d'r dh' ^ + E^,[n^] (74) 



47reo|r — -PI 

+ / rfV v{r)n\r) + E^,c[n^] . (75) 



i,'j&R 



The compensation charge density Z{f) = J2r ^r{^ is given as a sum of angular momentum 
dependent Gauss functions, which have an analytical plane-wave expansion. A similar term 
occurs also in the pseudopotential approach. In contrast to the norm-conserving pseudopotential 
approach, however, the compensation charge of an atom Zr is non-spherical and constantly 
adapts instantaneously to the environment. It is constructed such that 

n],{r) + Znir) - njj(r) - Z«(r) (76) 

has vanishing electrostatic multi-pole moments for each atomic site. With this choice, the elec- 
trostatic potentials of the augmentation densities vanish outside their spheres. This is the reason 
that there is no electrostatic interaction of the one-center parts between different sites. 
The compensation charge density as given here is still localized within the atomic regions. A 
technique similar to an Ewald summation, however, allows to replace it by a very extended 
charge density. Thus we can achieve, that the plane- wave convergence of the total energy is not 
affected by the auxiliary density. 

The potential v = J2r ^R^ which occurs in Eqs . 1731 and |75] enters the total energy in the form of 
"intelligent zeros" described in Eq. |7T] 

^ = Y1 [(i'nlvRli'n) - (^nl^fll^n)) = Y fn{i'n\vR\i'n) " ^ Eij{4>i\vR\(i>j) ■ (77) 
n n i,j&R 

The main reason for introducing this potential is to cancel the Coulomb singularity of the po- 
tential in the plane-wave part. The potential v allows to influence the plane-wave convergence 
beneficially, without changing the converged result, v must be localized within the augmenta- 
tion region, where Eq. l58l holds. 



Approximations 

Once the total energy functional provided in the previous section has been defined, everything 
else follows: Forces are partial derivatives with respect to atomic positions. The potential is the 
derivative of the non-kinetic energy contributions to the total energy with respect to the density. 
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and the auxiliary Hamiltonian follows from derivatives H\ipn) with respect to auxiliary wave 



functions. The fictitious-Lagrangian approach of Car and Parrinello ||98l does not allow any 
freedom in the way these derivatives are obtained. Anything else than analytic derivatives will 
violate energy conservation in a dynamical simulation. Since the expressions are straightfor- 
ward, even though rather involved, we will not discuss them here. 

All approximations are incorporated already in the total energy functional of the PAW method. 
What are those approximations? 

• Firstly we use the frozen-core approximation. In principle this approximation can be 



• The plane- wave expansion for the auxiliary wavefunctions must be complete. The plane- 
wave expansion is controlled easily by increasing the plane- wave cutoff defined as Epw = 
\h^G'^^^. Typically we use a plane-wave cutoff of 30 Ry. 

• The partial wave expansions must be converged. Typically we use one or two partial 
waves per angular momentum (£, m) and site. It should be noted that the partial wave 
expansion is not variational, because it changes the total energy functional and not the 
basis set for the auxiliary wavefunctions. 

Here, we do not discuss here numerical approximations such as the choice of the radial grid, 
since those are easily controlled. 

Relation to pseudopotentials 

We mentioned earlier that the pseudopotential approach can be derived as a well defined approx- 
imation from the PAW method: The augmentation part of the total energy = — for 
one atom is a functional of the one-center density matrix D defined in Eq. [68l The pseudopo- 
tential approach can be recovered if we truncate a Taylor expansion of ^E about the atomic 
density matrix after the linear term. The term linear in D is the energy related to the nonlocal 
pseudopotential. 



which can directly be compared to the total energy expression Eq. |49] of the pseudopotential 
method. The local potential v{r) of the pseudopotential approach is identical to the correspond- 
ing potential of the projector augmented-wave method. The remaining contributions in the PAW 
total energy, namely E, differ from the corresponding terms in Eq.|49]only in two features: our 
auxiliary density also contains an auxiliary core density, reflecting the nonlinear core correc- 
tion of the pseudopotential approach, and the compensation density Z{r) is non-spherical and 
depends on the wave function. Thus we can look at the PAW method also as a pseudopotential 



overcome. 




n 
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method with a pseudopotential that adapts instantaneously to the electronic environment. In the 
PAW method, the explicit nonlinear dependence of the total energy on the one-center density 
matrix is properly taken into account. 

What are the main advantages of the PAW method compared to the pseudopotential approach? 
Firstly all errors can be systematically controlled so that there are no transferability errors. As 
shown by Watson [|99l and Kresse HIOOH . most pseudopotentials fail for high spin atoms such as 
Cr. While it is probably true that pseudopotentials can be constructed that cope even with this 
situation, a failure can not be known beforehand, so that some empiricism remains in practice: 
A pseudopotential constructed from an isolated atom is not guaranteed to be accurate for a 
molecule. In contrast, the converged results of the PAW method do not depend on a reference 
system such as an isolated atom, because PAW uses the full density and potential. 
Like other all-electron methods, the PAW method provides access to the full charge and spin 
density, which is relevant, for example, for hyperfine parameters. Hyperfine parameters are 
sensitive probes of the electron density near the nucleus. In many situations they are the only 
information available that allows to deduce atomic structure and chemical environment of an 
atom from experiment. 

The plane-wave convergence is more rapid than in norm-conserving pseudopotentials and should 
in principle be equivalent to that of ultra-soft pseudopotentials lf88l . Compared to the ultra-soft 
pseudopotentials, however, the PAW method has the advantage that the total energy expression 
is less complex and can therefore be expected to be more efficient. 

The construction of pseudopotentials requires to determine a number of parameters. As they 
influence the results, their choice is critical. Also the PAW methods provides some flexibility 
in the choice of auxiliary partial waves. However, this choice does not influence the converged 
results. 

Recent developments 

Since the first implementation of the PAW method in the CP-PAW code ll58l . a number of 
groups have adopted the PAW method. The second implementation was done by the group 
of Holzwarth HlOlll . The resulting PWPAW code is freely available I1102II . This code is also 
used as a basis for the PAW implementation in the ABINIT project II104L An independent 
PAW code has been developed by Valiev and Weare II103L This implementation has entereed 
the NWChem code I1107II . An independent implementation of the PAW method is that of the 
VASP code niOOL The PAW method has also been implemented by W. Kromen II106II into the 
EStCoMPP code of Bliigel and Schroder. Another implementation is in the Quantum Espresso 
code II105II . A real-space-grid based version of the PAW method is the code GPAW developed 
by Mortensen et al. II108L 

Another branch of methods uses the reconstruction of the PAW method, without taking into 
account the full wavefunctions in the energy minimization. Following chemist notation, this 
approach could be termed "post-pseudopotential PAW". This development began with the eval- 
uation for hyperfine parameters from a pseudopotential calculation using the PAW reconstruc- 
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tion operator H 10911 and is now used in the pseudopotential approach to calculate properties that 
require the correct wavefunctions such as hyperfine parameters. 

The implementation of the PAW method by Kresse and Joubert HI 001 has been particularly use- 
ful as they had an implementation of PAW in the same code as the ultra-soft pseudopotentials, so 
that they could critically compare the two approaches with each other. Their conclusion is that 
both methods compare well in most cases, but they found that magnetic energies are seriously - 
by a factor two - in error in the pseudopotential approach, while the results of the PAW method 
were in line with other all-electron calculations using the linear augmented plane- wave method. 
As a short note, Kresse and Joubert incorrectly claim that their implementation is superior as 
it includes a term that is analogous to the non-linear core correction of pseudopotentials UllOII : 
this term however is already included in the original version in the form of the pseudized core 
density. 

Several extensions of the PAW have been done in the recent years: For applications in chemistry 
truly isolated systems are often of great interest. As any plane-wave based method introduces 
periodic images, the electrostatic interaction between these images can cause serious errors. 
The problem has been solved by mapping the charge density onto a point charge model, so that 
the electrostatic interaction could be subtracted out in a self-consistent manner [|111[| . In order 
to include the influence of the environment, the latter was simulated by simpler force fields 
using the molecular-mechanics-quantum-mechanics (QM-MM) approach II112L 

In order to overcome the limitations of the density functional theory several extensions have 
been performed. Bengone I1113II implemented the LDA-i-U approach into our CP-PAW code. 
Soon after this, Arnaud HI 1411 accomplished the implementation of the GW approximation into 
our CP-PAW code. The VASP-version of PAW [fTTSll and our CP-PAW code have now been 
extended to include a non-coUinear description of the magnetic moments. In a non-coUinear 
description, the Schrodinger equation is replaced by the Pauli equation with two-component 
spinor wavefunctions. 

The PAW method has proven useful to evaluate electric field gradients [I116II and magnetic hy- 
perfine parameters with high accuracy [I117II . Invaluable will be the prediction of NMR chemical 
shifts using the GIPAW method of Pickard and Mauri I1118L which is based on their earlier work 
111191 . While the GIPAW is implemented in a post-pseudopotential manner, the extension to a 
self-consistent PAW calculation should be straightforward. An post-pseudopotential approach 
has also been used to evaluate core level spectra HI 201 and momentum matrix elements [|121| . 
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Appendices 

A Model exchange-correlation energy 

We consider a model with a constant density, and a hole function, which describes a situation, 
where all electrons of the same spin are repelled completely from a sphere centered at the 
reference electron 
The hole function has the form 



/i(f, fo) 



\n{ro) for |r — ro| < rt 



2 

otherwise 



where n{r) is the electron density and the hole radius = y is the radius of the sphere 
which is determined such that the exchange correlation hole integrates to —1, i.e. ^rli^^n) = 
1. 

The potential of a homogeneously charged sphere with radius r/j and one positive charge is 



v{r) 



4vreo 1 _i 



r 



for r > Th 



where r = \ f — fo\. 

With Eq.[T2lwe obtain for the potential contribution of the exchange correlation energy 



3 J2n 



4 



Uxc = - / d'^r n(r)v(r = 0) = - d^r \ 

J ^ ' ^ ' J 47reo4V 3 

B Large-gradient limit of the enhancement factor 

An exponentially decaying density 

n(r) = exp(— Ar) (79) 

has a reduced gradient 

I Vnl 1 
a; := = Aexp(+-Ar) (80) 

77,3 3 

We make the following ansatz for the exchange correlation energy per electron 

exc{n,x) = -Cn^F{x) (81) 
where only the local exchange has been used and C is a constant. 

Enforcing the long-distance limit of the exchange correlation energy per electron for exponen- 
tially decaying densities 

e.,((n(r),x(r)) = -l^ (82) 
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yields 



F{x) = , 

47reor(x)2C?i3 (r(x)) 

Using Eqs.|79]and[80l we express the radius and the density by the reduced gradient, i.e. 

r{x) 
n{x) 



--(^ln[A] - ln[x] 
n(r(x)) = A^x~^ 



(83) 

(84) 
(85) 



and obtain 

Fix) 



47ren 



x—^oo 



-f( ln[A] -ln[x] 

x'^ 



2CXx' 



( 







X 



Aneo ■ 6C x ln(A) — x \n{x) 



(86) 



^47reo • 6C / x\n{x) 

Now we need to ensure that -F(O) = 1 so that the gradient correction vanishes for the homo- 
geneous electron gas, and that F(x) = F{—x) to enforce spin reversal symmetry. There are 
several possible interpolations for these requirements, but the most simple one is 

Px'^ 



F(x) = 1 

^ ^ 1 + ^ ■ 6C/3a; ■ asinh(x) 

This is the enhancement factor for exchange used by Becke in his B88 functional |[24l 



(87) 
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